// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//  File: avtTriangulateRegularPointsFilter.C
// ************************************************************************* //

#include <avtTriangulateRegularPointsFilter.h>

#include <vtkCellArray.h>
#include <vtkFloatArray.h>
#include <vtkPolyData.h>

#include <DebugStream.h>

// ****************************************************************************
//  Method: avtTriangulateRegularPointsFilter constructor
//
//  Programmer: brugger1 -- generated by xml2avt
//  Creation:   Fri Jun 5 17:15:09 PST 2009
//
// ****************************************************************************

avtTriangulateRegularPointsFilter::avtTriangulateRegularPointsFilter()
{
}


// ****************************************************************************
//  Method: avtTriangulateRegularPointsFilter destructor
//
//  Programmer: brugger1 -- generated by xml2avt
//  Creation:   Fri Jun 5 17:15:09 PST 2009
//
//  Modifications:
//
// ****************************************************************************

avtTriangulateRegularPointsFilter::~avtTriangulateRegularPointsFilter()
{
}


// ****************************************************************************
//  Method:  avtTriangulateRegularPointsFilter::Create
//
//  Programmer: brugger1 -- generated by xml2avt
//  Creation:   Fri Jun 5 17:15:09 PST 2009
//
// ****************************************************************************

avtFilter *
avtTriangulateRegularPointsFilter::Create()
{
    return new avtTriangulateRegularPointsFilter();
}


// ****************************************************************************
//  Method:      avtTriangulateRegularPointsFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: brugger1 -- generated by xml2avt
//  Creation:   Fri Jun 5 17:15:09 PST 2009
//
// ****************************************************************************

void
avtTriangulateRegularPointsFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const TriangulateRegularPointsAttributes*)a;
}


// ****************************************************************************
//  Method: avtTriangulateRegularPointsFilter::Equivalent
//
//  Purpose:
//      Returns true if creating a new avtTriangulateRegularPointsFilter with the given
//      parameters would result in an equivalent avtTriangulateRegularPointsFilter.
//
//  Programmer: brugger1 -- generated by xml2avt
//  Creation:   Fri Jun 5 17:15:09 PST 2009
//
// ****************************************************************************

bool
avtTriangulateRegularPointsFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(TriangulateRegularPointsAttributes*)a);
}


// ****************************************************************************
//  Method: avtTriangulateRegularPointsFilter::ExecuteData
//
//  Purpose:
//      Sends the specified input and output through the TriangulateRegularPoints filter.
//
//  Arguments:
//      in_dr      The input data representation.
//
//  Returns:       The output data representation.
//
//  Programmer: brugger1 -- generated by xml2avt
//  Creation:   Fri Jun 5 17:15:09 PST 2009
//
//  Modifications:
//    Eric Brugger, Wed Jun 17 17:19:40 PDT 2009
//    Replaced several calls to the round function with appropriate
//    substitutes.
//
//    Eric Brugger, Tue Aug 19 09:41:40 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
// ****************************************************************************

avtDataRepresentation *
avtTriangulateRegularPointsFilter::ExecuteData(avtDataRepresentation *in_dr)
{
    //
    // Get the VTK data set.
    //
    vtkDataSet *in_ds = in_dr->GetDataVTK();

    //
    // The input should be a point mesh, so its type is poly data.
    //
    if (in_ds->GetDataObjectType() != VTK_POLY_DATA)
    {
        EXCEPTION1(ImproperUseException,
                   "Expected a vtkPolyData in the avtTriangulateRegularPointsFilter.");
    }
    vtkPolyData *in_pd = vtkPolyData::SafeDownCast(in_ds);

    int nPoints = in_pd->GetNumberOfPoints();
    float *points = vtkFloatArray::SafeDownCast(in_pd->GetPoints()->GetData())->GetPointer(0);
    if (points == NULL)
    {
        EXCEPTION1(ImproperUseException, "Expected float coordinate data.");
    }

    //
    // Get the attributes.
    //
    bool useXSpacing = atts.GetUseXGridSpacing();
    bool useYSpacing = atts.GetUseYGridSpacing();
    double xSpacing = atts.GetXGridSpacing();
    double ySpacing = atts.GetYGridSpacing();

    //
    // Calculate x and y extents.
    //
    double xmin = points[0];
    double xmax = points[0];
    double ymin = points[1];
    double ymax = points[1];
    for (int j = 1; j < nPoints; j++)
    {
        xmin = (points[j*3+0] < xmin) ? points[j*3+0] : xmin;
        xmax = (points[j*3+0] > xmax) ? points[j*3+0] : xmax;
        ymin = (points[j*3+1] < ymin) ? points[j*3+1] : ymin;
        ymax = (points[j*3+1] > ymax) ? points[j*3+1] : ymax;
    }
    debug1 << "xmin=" << xmin << ",xmax=" << xmax << endl;
    debug1 << "ymin=" << ymin << ",ymax=" << ymax << endl;

    //
    // Determine the grid spacing and dimensions.
    //
    double dx = (xmax - xmin);
    double dy = (ymax - ymin);
    double dxMin = dx / 2000.;
    double dyMin = dy / 2000.;
    for (int j = 0; j < nPoints; j++)
    {
        double dxTmp = points[j*3+0] - xmin;
        double dyTmp = points[j*3+1] - ymin;
        dx = (dxTmp < dx && dxTmp > dxMin) ? dxTmp : dx;
        dy = (dyTmp < dy && dyTmp > dyMin) ? dyTmp : dy;
    }
    debug1 << "dx=" << dx << ",dy=" << dy << endl;

    if (useXSpacing) dx = xSpacing;
    if (useYSpacing) dy = ySpacing;

    //
    // Calculate the grid dimensions.  The 0.5 is being added to round the
    // values.  The round function was not used since it doesn't exist on
    // Windows.
    //
    int nx = (int) ((xmax - xmin) / dx + 1. + 0.5);
    int ny = (int) ((ymax - ymin) / dy + 1. + 0.5);

    //
    // Map each point to the grid. The 0.5 is being added to ix and iy
    // in the for loop to round the values.
    //
    int *rgrid = new int[nx*ny];
    for (int i = 0; i < nx*ny; i++)
        rgrid[i] = -1;

    int imax = -1;
    for (int i = 0 ; i < nPoints; i++)
    {
        int ix = (int) ((points[i*3+0] - xmin) / dx + 0.5);
        int iy = (int) ((points[i*3+1] - ymin) / dy + 0.5);
        if (iy*nx+ix > imax) imax = iy*nx+ix;
        rgrid[iy*nx+ix] = i;
    }

    //
    // Create the triangles from the grid.
    //
    vtkIdType ptIds[3];
    vtkCellArray *triangles = vtkCellArray::New();
    for (int i = 0; i < nx-1; i++)
    {
        for (int j = 0; j < ny-1; j++)
        {
            int ndx = ((rgrid[j*nx+i] > 0) ? 1 : 0) +
                      ((rgrid[j*nx+i+1] > 0) ? 2 : 0) +
                      ((rgrid[(j+1)*nx+i] > 0) ? 4 : 0) +
                      ((rgrid[(j+1)*nx+i+1] > 0) ? 8 : 0);
            switch (ndx)
            {
              case 7:
                ptIds[0] = rgrid[j*nx+i];
                ptIds[1] = rgrid[(j+1)*nx+i];
                ptIds[2] = rgrid[j*nx+i+1];
                triangles->InsertNextCell(3, ptIds);
                break;
              case 11:
                ptIds[0] = rgrid[(j+1)*nx+i+1];
                ptIds[1] = rgrid[j*nx+i];
                ptIds[2] = rgrid[j*nx+i+1];
                triangles->InsertNextCell(3, ptIds);
                break;
              case 13:
                ptIds[0] = rgrid[(j+1)*nx+i];
                ptIds[1] = rgrid[(j+1)*nx+i+1];
                ptIds[2] = rgrid[j*nx+i];
                triangles->InsertNextCell(3, ptIds);
                break;
              case 14:
                ptIds[0] = rgrid[(j+1)*nx+i];
                ptIds[1] = rgrid[j*nx+i+1];
                ptIds[2] = rgrid[(j+1)*nx+i+1];
                triangles->InsertNextCell(3, ptIds);
                break;
              case 15:
                ptIds[0] = rgrid[j*nx+i];
                ptIds[1] = rgrid[(j+1)*nx+i];
                ptIds[2] = rgrid[j*nx+i+1];
                triangles->InsertNextCell(3, ptIds);
                ptIds[0] = rgrid[(j+1)*nx+i];
                ptIds[1] = rgrid[j*nx+i+1];
                ptIds[2] = rgrid[(j+1)*nx+i+1];
                triangles->InsertNextCell(3, ptIds);
                break;
              default:
                break;
            }
        }
    }

    delete [] rgrid;

    //
    // Create the output data set.
    //
    vtkPolyData *out_pd = in_pd->NewInstance();
    out_pd->ShallowCopy(in_pd);

    out_pd->SetVerts(NULL);
    out_pd->SetPolys(triangles);
    triangles->Delete();

    //
    // Return the output data set.
    //
    avtDataRepresentation *out_dr = new avtDataRepresentation(out_pd,
        in_dr->GetDomain(), in_dr->GetLabel());

    out_pd->Delete();

    return out_dr;
}
